
#include "CFDProblem.h"
#include "Eigen/Geometry"
#include "NonlinearSystem.h"
using Eigen::Vector3d;

template<>
InputParameters validParams<CFDProblem>()
{
  InputParameters params = validParams<FEProblem>();
  MooseEnum vis_type(CFDProblem::viscousType());
  params.addParam<MooseEnum>("vis_type", vis_type, "粘性计算方法");
  MooseEnum flux_type(CFDProblem::fluxRiemannType());
  params.addParam<MooseEnum>("flux_type", flux_type, "Riemann通量类型");
  params.addParam<Real>("mach",  0.1, "马赫数");
  params.addParam<Real>("gamma", 1.4, "比热比");
  params.addParam<Real>("reynolds", 1, "雷诺数");
  params.addParam<Real>("prandtl", 0.72, "prandtl数");
  params.addParam<Real>("attack", 0., "攻角");
  params.addParam<Real>("sideslip", 0., "侧滑角");
  params.addParam<Real>("pitch", 0., "俯仰角");
  params.addParam<Real>("yaw", 180., "偏航角");
  params.addParam<Real>("roll", -90., "滚转角");
  params.addParam<Real>("init_vel", 1.0, "初始流体速度");
  params.addParam<unsigned int>("jacobian_delay", 1, "jacobian矩阵更新频率");
  return params;
}

CFDProblem::CFDProblem(const InputParameters &params) :
	FEProblem(params),
//	_var_order(1),
	_vis_type(getParam<MooseEnum>("vis_type")),
	_flux_type(getParam<MooseEnum>("flux_type")),
	_mach(getParam<Real>("mach")),
	_gamma(getParam<Real>("gamma")),
	_reynolds(getParam<Real>("reynolds")),
	_prandtl(getParam<Real>("prandtl")),

	_attack(getParam<Real>("attack")*libMesh::pi/180),
	_sideslip(getParam<Real>("sideslip")*libMesh::pi/180),
	_pitch(getParam<Real>("pitch")*libMesh::pi/180),
	_yaw(getParam<Real>("yaw")*libMesh::pi/180),
	_roll(getParam<Real>("roll")*libMesh::pi/180),

	_attitude(_attack, _sideslip, _pitch, _yaw, _roll),
    _velocity(getParam<Real>("init_vel")),
	_jacobian_delay(getParam<unsigned int>("jacobian_delay"))
{
}

void CFDProblem::valueExact(Real* value, Real t, const Point& p)
{
	mooseError("CFDProblem::valueExact 不可调用，在子类实现");
}

void CFDProblem::initialSetup()
{
	FEProblem::initialSetup();
	computeIndicatorsAndMarkers();
}

void CFDProblem::computeJacobian(const NumericVector<Number> & soln, SparseMatrix<Number> &  jacobian)
{
	if(_t_step == 1|| _t_step == 2 || (_t_step-1) % _jacobian_delay == 0 || _app.isRecovering() || _app.isRestarting() || !converged())
		Parent::computeJacobian(soln, jacobian);
}

MooseEnum CFDProblem::viscousType()
{
  return MooseEnum("INVISCOUS CONSTANT  SUTHERLAND", "CONSTANT");
}

MooseEnum CFDProblem::fluxRiemannType()
{
//  return MooseEnum("Lax-F HLL  HLLC-PV HLLC-Roe", "Lax-F");
  return MooseEnum("Lax-F HLL", "Lax-F");
}

void CFDProblem::primaryTOConservation(Real *pri, Real * cons)
{
	Real rho = pri[0];
	Real u = pri[1];
	Real v = pri[2];
	Real w = pri[3];
	Real p = pri[4];

	cons[0] = rho;
	cons[1] = rho*u;
	cons[2] = rho*v;
	cons[3] = rho*w;
	cons[4] = p/(_gamma-1) +0.5*rho*(u*u+v*v+w*w);
}

void CFDProblem::initialCondition(Real *value, Real t, const Point &p)
{
//	Real density = 1;
//	Vector3d vel  = _velocity*(_attitude.earthFromWind()*Vector3d::UnitX());
//	if(_mesh.dimension() < 3)
//		vel(2) = 0;
//	if(_mesh.dimension() < 2)
//		vel(1) = 0;
//	Vector3d mom = density*vel;
//
//	Real pre = 1./_gamma/_mach/_mach;
//	Real rhoe = pre/(_gamma-1) + 0.5*density*(_velocity*_velocity);
//	value[0] = density;
//	value[1] = mom(0);
//	value[2] = mom(1);
//	value[3] = mom(2);
//	value[4] = rhoe;
	mooseError("CFDProblem::initialCondition 不可调用，改用CFDPassFlowIC");
	valueExact(value, t, p);
}

void CFDProblem::boundaryCondition(Real *value, Real t, const Point &p)
{
	valueExact(value, t, p);
}
